
# 5 新发感染区县增加的情况------

yearly_incidence_county <- readxl::read_excel("./output/表格3.1.1_发病数和发病率_逐年.xlsx")

# 5.1.1 2005-2023年历年新发感染区县------
unique(yearly_incidence_county$city_code)
# map_data$city_code %>% unique()

yearly_incidence_county$county_code %>% unique() %in% "654223" #650502

# yearly_incidence_county %>% filter(city_code == "6542") %>% View()
# 获取所有年份和区县的组合，并初始化状态列
full_grid <- expand.grid(year = seq(min(yearly_incidence_county$year), 
                                    max(yearly_incidence_county$year)),
                         county_code = unique(yearly_incidence_county$county_code)) %>%
  left_join(yearly_incidence_county, by = c("year", "county_code")) %>%
  mutate(status = "无病例报道")


# 更新状态列
update_status <- function(df) {
  df %>%
    group_by(county_code) %>%
    arrange(year) %>%
    mutate(
      cumulative_infections = cumsum(!is.na(cases) & cases > 0),
      status = case_when(
        !is.na(cases) & cases > 0 & cumulative_infections == 1 ~ "当年首次感染",
        !is.na(cases) & cases > 0 ~ "曾经感染过",
        is.na(cases) | cases == 0 ~ "无病例报道"
      )
    ) %>%
    ungroup() %>%
    mutate(status = ifelse(is.na(cumulative_infections), "从来没有感染过", status)) %>%
    select(-cumulative_infections)
}

full_grid <- update_status(full_grid)
full_grid %>% filter(year == 2023)
writexl::write_xlsx(full_grid, "./output/表格5_emerging_infection_status.xlsx")
full_grid$status %>% unique()

# 3.3.2 统计逐年状态数 -------
# 计算每年的状态数量
status_counts <- full_grid %>%
  group_by(year, status) %>%
  summarise(count = n(), .groups = 'drop') %>%
  complete(year, status, fill = list(count = 0))
summary(status_counts)

status_counts <- status_counts %>%
  group_by(year) %>%
  arrange(desc(status)) %>% # 按照新设定的因子顺序排序
  mutate(cumulative_count = cumsum(count),
         label_position = lag(cumulative_count, default = 0) + count / 2) # 标签位于中间
writexl::write_xlsx(status_counts, "./output/表格5_status_counts.xlsx")

# if (!"感染状态_计数" %in% names(wb2)) {addWorksheet(wb2, "感染状态_计数")}
# writeData(wb2, sheet = "感染状态_计数", x = status_counts)

# 绘制堆叠柱状图并添加所有状态的标签
ggplot(status_counts, aes(x = factor(year), y = count, fill = status)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(y = label_position, label = count),
            vjust = 0.5, color = "black", size = 3, hjust = 0.5) + # 调整vjust使标签居中
  scale_fill_manual(values = c("从来没有感染过" = "lightgreen",
                               "曾经感染过" = "orange",
                               "当年首次感染" = "red",
                               "无病例报道" = "gray80"),
                    na.value = "gray70") +
  labs(title = "Infection Status of Counties Over Years",
       x = "Year",
       y = "Number of Counties",
       fill = "Status") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) # 旋转x轴标签以提高可读性
# ggsave("./figs/图3.3.2 逐年状态区县数分布.png")

# 5.1.2 绘图_逐年------


library(dplyr)
library(ggplot2)
library(scales)

# 步骤1: 计算年度报告区县数和比例
total_counties <- n_distinct(yearly_incidence_county$county_code)

yearly_report_ratio <- yearly_incidence_county %>%
  group_by(year) %>%
  summarise(
    reporting_counties = sum(cases > 0, na.rm = TRUE),
    reporting_ratio = reporting_counties / total_counties
  )

# 步骤2: 计算累计受影响区县数
cumulative_counties <- yearly_incidence_county %>%
  filter(cases > 0) %>%
  group_by(county_code) %>%
  summarise(first_year = min(year)) %>%
  group_by(first_year) %>%
  summarise(new_counties = n()) %>%
  complete(first_year = 2005:2023, fill = list(new_counties = 0)) %>%
  arrange(first_year) %>%
  mutate(cumulative_counties = cumsum(new_counties)) %>%
  rename(year = first_year)

# 步骤3: 合并数据
combined_data <- yearly_report_ratio %>%
  left_join(cumulative_counties, by = "year")

# 步骤4: 绘制极简风格双轴图
ggplot(combined_data, aes(x = year)) +
  # 柱状图（累计区县数）
  geom_col(
    aes(y = cumulative_counties, fill = "Cumulative Affected Counties"),
    width = 0.7,
    alpha = 0.9,
    show.legend = TRUE
  ) +
  # 柱顶标签（年度报告数）
  geom_text(
    aes(y = cumulative_counties, label = reporting_counties),
    vjust = -0.7,
    color =  "orange",# "#1f77b4", 
    size = 3.8,
    fontface = "bold"
  ) +
  # 折线图（年度报告比例）
  geom_line(
    aes(
      y = reporting_ratio * max(cumulative_counties),
      color = "Annual Reporting Proportion"
    ),
    linewidth = 1.1,
    linetype = "longdash"
  ) +
  geom_point(
    aes(
      y = reporting_ratio * max(cumulative_counties),
      color = "Annual Reporting Proportion"
    ),
    size = 3.2,
    shape = 18
  ) +
  # 双轴设置
  scale_y_continuous(
    name = "Cumulative Affected Counties",
    sec.axis = sec_axis(
      ~ . / max(combined_data$cumulative_counties),
      name = "County Reporting Proportion",
      labels = percent_format(accuracy = 1)
    ),
    expand = expansion(mult = c(0, 0.12))
  ) +
  # 颜色与图例
  scale_fill_manual(
    name = "",
    values = c("Cumulative Affected Counties" = "orange")#"#1f77b4"
  ) +
  scale_color_manual(
    name = "",
    values = c("Annual Reporting Proportion" = "#d62728")
  ) +
  # 极简主题
  theme_classic() +
  theme(
    axis.title.y.left = element_text(
      color = "orange",#"#1f77b4"
      size = 12,
      margin = margin(r = 8)
    ),
    axis.title.y.right = element_text(
      color = "#d62728",
      size = 12,
      margin = margin(l = 8)
    ),
    axis.text.x = element_text(
      angle = 45,
      hjust = 1,
      vjust = 1.1,
      size = 10,
      color = "black"
    ),
    axis.line = element_line(color = "black", linewidth = 0.4),
    axis.ticks = element_line(color = "black", linewidth = 0.4),
    legend.position = "top",
    legend.direction = "horizontal",
    legend.spacing.x = unit(0.4, "cm"),
    legend.text = element_text(size = 11),
    panel.background = element_rect(fill = "white"),
    plot.margin = margin(10, 15, 10, 15)
  ) +
  guides(
    fill = guide_legend(order = 1, title.position = "top"),
    color = guide_legend(order = 2, title.position = "top")
  ) +
  labs(x = NULL)

# 保存高清图（调整dpi和尺寸）
ggsave("./主图/Fig_感染区县柱状图.png", width = 8, height = 6, dpi = 600)

# library(dplyr)
# library(ggplot2)
# library(scales)
# 
# # 步骤1: 计算年度报告区县数和比例
# total_counties <- n_distinct(yearly_incidence_county$county_code)
# 
# yearly_report_ratio <- yearly_incidence_county %>%
#   group_by(year) %>%
#   summarise(
#     reporting_counties = sum(cases > 0, na.rm = TRUE),  # 年度报告区县数
#     reporting_ratio = reporting_counties / total_counties
#   )
# 
# # 步骤2: 计算累计受影响区县数
# cumulative_counties <- yearly_incidence_county %>%
#   filter(cases > 0) %>%
#   group_by(county_code) %>%
#   summarise(first_year = min(year)) %>%
#   group_by(first_year) %>%
#   summarise(new_counties = n()) %>%
#   complete(first_year = 2005:2023, fill = list(new_counties = 0)) %>%
#   arrange(first_year) %>%
#   mutate(cumulative_counties = cumsum(new_counties)) %>%
#   rename(year = first_year)
# 
# # 步骤3: 合并数据（保留reporting_counties）
# combined_data <- yearly_report_ratio %>%
#   left_join(cumulative_counties, by = "year")
# 
# # 步骤4: 绘制双轴图（柱状图 + 数值标签 + 折线）
# ggplot(combined_data, aes(x = year)) +
#   # 累计区县数（柱状图）
#   geom_col(
#     aes(y = cumulative_counties, fill = "Cumulative Affected Counties"),
#     width = 0.7,
#     alpha = 0.8
#   ) +
#   # 标注年度报告区县数（柱顶）
#   geom_text(
#     aes(
#       y = cumulative_counties,
#       label = reporting_counties  # 直接显示reporting_counties值
#     ),
#     vjust = -0.5,                # 垂直位置（柱顶上方）
#     color = "#1f77b4",
#     size = 3.5,
#     fontface = "bold"
#   ) +
#   # 年度报告率（折线）
#   geom_line(
#     aes(
#       y = reporting_ratio * max(cumulative_counties),
#       color = "Annual Reporting Rate"
#     ),
#     linewidth = 1.2,
#     linetype = "dashed"
#   ) +
#   geom_point(
#     aes(
#       y = reporting_ratio * max(cumulative_counties),
#       color = "Annual Reporting Rate"
#     ),
#     size = 3
#   ) +
#   # 双轴设置
#   scale_y_continuous(
#     name = "Cumulative Affected Counties",
#     sec.axis = sec_axis(
#       ~ . / max(combined_data$cumulative_counties),
#       name = "County Reporting Proportion",
#       labels = percent_format(accuracy = 1)
#     ),
#     expand = expansion(mult = c(0, 0.1))  # 扩展顶部空间以适应标签
#   ) +
#   # 颜色与图例
#   scale_fill_manual(
#     name = "",
#     values = c("Cumulative Affected Counties" = "#1f77b4")
#   ) +
#   scale_color_manual(
#     name = "",
#     values = c("Annual Reporting Rate" = "#d62728")
#   ) +
#   # 主题优化
#   labs(x = NULL) +#"Year"
#   theme_minimal() +
#   theme(
#     axis.title.y.left = element_text(
#       color = "#1f77b4",
#       size = 12,
#       margin = margin(r = 10)
#     ),
#     axis.title.y.right = element_text(
#       color = "#d62728",
#       size = 12,
#       margin = margin(l = 10)
#     ),
#     axis.text.x = element_text(angle = 45, hjust = 1, size = 10),
#     legend.position = "top",
#     legend.spacing.x = unit(0.3, "cm"),
#     panel.grid.major = element_line(color = "grey90"),
#     panel.grid.minor = element_blank()
#   ) +
#   guides(
#     fill = guide_legend(order = 1),
#     color = guide_legend(order = 2)
#   )
# # 保存高清图
# ggsave("./主图/Fig 感染区县柱状图.png",dpi = 600)

# 5.1.3 Kimi的关于统计“首次报告"的计算算法：用min函数统计每个区县的min，然后转换数据

library(dplyr)
library(ggplot2)
library(scales)

# 1. 逐年的有布病病例报告的区县数目
total_counties <- n_distinct(yearly_incidence_county$county_code)  # 总区县数

yearly_reported_counties <- yearly_incidence_county %>%
  group_by(year) %>%
  summarise(
    reported_counties = sum(cases > 0, na.rm = TRUE)
  ) %>%
  mutate(
    reported_ratio = reported_counties / total_counties
  )

# 2. 自2005年以来，每年的新出现有病例感染的区县数目
first_report_year <- yearly_incidence_county %>%
  filter(cases > 0) %>%
  group_by(county_code) %>%
  summarise(first_report_year = min(year)) %>%
  ungroup()

new_counties_per_year <- first_report_year %>%
  count(first_report_year) %>%
  rename(year = first_report_year) %>%
  complete(year = 2005:2023, fill = list(n = 0)) %>%
  arrange(year) %>%
  mutate(cumulative_affected = cumsum(n))

# 3. 合并数据
combined_data <- yearly_reported_counties %>%
  left_join(new_counties_per_year, by = "year")

# # 4. 绘制双轴图
# ggplot(combined_data, aes(x = year)) +
#   # 累计受影响区县数（左轴）- 柱状图
#   geom_col(aes(y = cumulative_affected, fill = "Cumulative Affected Counties"), 
#            width = 0.6, alpha = 0.7) +
#   # 报告比例（右轴）
#   geom_line(aes(y = reported_ratio * max(cumulative_affected), 
#                 color = "County Reporting Proportion"), 
#             linewidth = 1.2, linetype = "dashed") +
#   geom_point(aes(y = reported_ratio * max(cumulative_affected), 
#                  color = "County Reporting Proportion"), 
#              size = 2.5) +
#   # 双轴设置
#   scale_y_continuous(
#     name = "Cumulative Affected Counties",
#     sec.axis = sec_axis(
#       ~ . / max(combined_data$cumulative_affected), 
#       name = "County Reporting Proportion",
#       labels = percent_format(accuracy = 1)
#     ),
#     limits = c(0, max(combined_data$cumulative_affected) * 1.05)
#   ) +
#   # 颜色与图例
#   scale_color_manual(
#     name = "Metrics",
#     values = c(
#       "Cumulative Affected Counties" = "#1f77b4",
#       "County Reporting Proportion" = "#d62728"
#     ),
#     guide = guide_legend(override.aes = list(linetype = c("blank", "dashed")))
#   ) +
#   scale_fill_manual(
#     name = "Metrics",
#     values = c(
#       "Cumulative Affected Counties" = "#1f77b4"
#     )
#   ) +
#   # 标签与主题
#   labs(
#     x = "Year"
#   ) +
#   theme_minimal() +
#   theme(
#     axis.title.y.left = element_text(color = "#1f77b4", size = 12),
#     axis.title.y.right = element_text(color = "#d62728", size = 12),
#     axis.text.x = element_text(angle = 45, hjust = 1),
#     legend.position = "bottom",
#     legend.title = element_blank(),
#     panel.grid.major = element_line(color = "grey90"),
#     panel.grid.minor = element_blank()
#   )

# 5.2 分阶段新发感染区县------
# 
# 定义时间段
time_periods <- list(
  all_years = 2004:2023,
  y2005_2015 = 2005:2015,
  y2016_2020 = 2016:2020,
  y2021_2023 = 2021:2023
)

library(dplyr)
library(ggplot2)
library(sf)
library(writexl)
# library(ggsave)
library(stringr)

# 初始化函数以更新状态列
update_status <- function(df) {
  df %>%
    group_by(county_code) %>%
    arrange(year) %>%
    mutate(
      cumulative_infections = cumsum(!is.na(cases) & cases > 0),
      status = case_when(
        !is.na(cases) & cases > 0 & cumulative_infections == 1 ~ "当年首次感染",
        !is.na(cases) & cases > 0 ~ "曾经感染过",
        TRUE ~ "无病例报道" # 使用TRUE来捕获所有其他情况
      ),
      status = ifelse(is.na(cumulative_infections), "从来没有感染过", status)
    ) %>%
    ungroup() %>%
    select(-cumulative_infections)
}

# 获取所有年份和区县的组合，并初始化状态列
full_grid <- expand.grid(
  year = seq(min(yearly_incidence_county$year), max(yearly_incidence_county$year)),
  county_code = unique(yearly_incidence_county$county_code)
) %>%
  left_join(yearly_incidence_county, by = c("year", "county_code")) %>%
  mutate(status = "无病例报道") %>%
  update_status()

# 保存 full_grid 数据
write_xlsx(full_grid, "./output/表格5.1_感染状态_历年区县.xlsx")

# 绘制逐年地图函数
plot_status_map <- function(rate_data, title, map_data) {
  rate_map_data <- map_data %>%
    left_join(rate_data, by = c("GB1999" = "county_code"))
  
  centroids <- st_centroid(st_geometry(rate_map_data))
  rate_map_data$centroid_x <- st_coordinates(centroids)[,1]
  rate_map_data$centroid_y <- st_coordinates(centroids)[,2]
  
  ggplot() +
    geom_sf(data = rate_map_data, aes(fill = status)) +
    scale_fill_manual(values = c("从来不曾感染过" = "lightgreen",
                                 "曾经感染过" = "orange",
                                 "当年首次感染" = "red",
                                 "无病例报道" = "white"),
                      na.value = "gray50") +
    labs(title = title, fill = "Status") +
    theme_void() + 
    theme(
      plot.background = element_rect(fill = "white"),  
      plot.title = element_text(hjust = 0.5),          
      legend.position = "bottom",                      
      legend.box.background = element_rect(fill = "white", color = NA)  
    ) +
    annotation_north_arrow(location = "tl", which_north = "true",
                           style = north_arrow_fancy_orienteering) + 
    annotation_scale(location = "bl") + 
    geom_text(data = rate_map_data, aes(x = centroid_x, y = centroid_y, label = NAME), size = 3, vjust = -0.5)
}

# 分别绘制不同年份的地图并保存图片
dir.create("figs/逐年状态地图", showWarnings = FALSE, recursive = TRUE)
plots <- lapply(unique(full_grid$year), function(yr) {
  year_data <- filter(full_grid, year == yr)
  
  p <- plot_status_map(year_data, title = paste("Infection Status for", yr), map_data = map_data)
  
  ##ggsave(filename = file.path("figs/逐年状态地图", paste0("infection_status_", yr, ".png")),
  # plot = p, width = 10, height = 8, dpi = 300)
  
  p
})

# 计算每年的状态数量
status_counts <- full_grid %>%
  count(year, status, name = "count") %>%
  complete(year, status, fill = list(count = 0)) %>%
  group_by(year) %>%
  arrange(desc(status)) %>%
  mutate(cumulative_count = cumsum(count),
         label_position = lag(cumulative_count, default = 0) + count / 2)

# 保存每年状态数量
# status_counts %>% View()
write_xlsx(status_counts, "./output/表格5.1_感染状态_新发感染区县数量.xlsx")

# 定义所有可能的状态
status_counts$status %>% unique()
status_levels <- c("从来不曾感染过", "无病例报道", "曾经感染过", "当年首次感染")

# 绘制逐年堆叠柱状图
ggplot(status_counts, aes(x = factor(year), y = count, fill = status)) +
  geom_bar(stat = "identity", position = "stack") +
  geom_text(aes(y = label_position, label = count), vjust = 0.5, color = "black", size = 3, hjust = 0.5) +
  scale_fill_manual(values = c("从来不曾感染过" = "lightgreen",
                               "曾经感染过" = "orange",
                               "当年首次感染" = "red",
                               "无病例报道" = "gray90"),
                    na.value = "gray50",
                    limits = status_levels) +
  labs(title = "Infection Status of Counties Over Years",
       x = "Year",
       y = "Number of Counties",
       fill = "Status") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
ggsave("./figs/图3.3.2 逐年状态区县数分布.png")

# 分析并绘制各时间段的地图
dir.create("figs/分阶段状态地图", showWarnings = FALSE, recursive = TRUE)

plots_by_period <- lapply(names(time_periods), function(period_name) {
  years_in_period <- time_periods[[period_name]]
  
  period_data <- full_grid %>%
    filter(year %in% years_in_period) %>%
    update_status() %>%
    distinct(county_code, .keep_all = TRUE)
  
  p <- plot_status_map(period_data, title = paste("Infection Status for", period_name), map_data = map_data)
  
  ##ggsave(filename = file.path("figs/分阶段状态地图", paste0("infection_status_", period_name, ".png")),
  # plot = p, width = 10, height = 8, dpi = 300)
  p
})

# 每个阶段的状态数量
status_counts_combined <- bind_rows(lapply(names(time_periods), function(period_name) {
  years_in_period <- time_periods[[period_name]]
  
  period_data <- full_grid %>%
    filter(year %in% years_in_period) %>%
    update_status() %>%
    distinct(county_code, .keep_all = TRUE) %>%
    count(status, name = "count") %>%
    complete(status, fill = list(count = 0)) %>%
    mutate(period = period_name)
}))

# 保存结果
status_counts_combined
write_xlsx(status_counts_combined, "./output/表格5_感染状态_分阶段.xlsx")


# if (!"感染状态_计数_分阶段" %in% names(wb2)) {addWorksheet(wb2, "感染状态_计数_分阶段")}
# writeData(wb2, sheet = "感染状态_计数_分阶段", x = status_counts_combined)

# saveWorkbook(wb1, file = "./output/0_发病率.xlsx", overwrite = TRUE)
# saveWorkbook(wb2, file = "./output/0_莫兰指数.xlsx", overwrite = TRUE)

# # 绘制分阶段堆叠柱状图
# ggplot(status_counts_combined, aes(x = factor(period), y = count, fill = status)) +
#   geom_bar(stat = "identity", position = "stack") +
#   geom_text(aes(y = count / 2, label = count), vjust = 0.5, color = "black", size = 3, hjust = 0.5) +
#   scale_fill_manual(values = c("从来不曾感染过" = "lightgreen",
#                                "曾经感染过" = "orange",
#                                "当年首次感染" = "red",
#                                "无病例报道" = "gray80"),
#                     na.value = "gray50",
#                     limits = status_levels) +
#   labs(title = "Infection Status of Counties Over Time Periods",
#        x = "Time Period",
#        y = "Number of Counties",
#        fill = "Status") +
#   theme_minimal() +
#   theme(axis.text.x = element_text(angle = 45, hjust = 1)) 
# 
##ggsave("./figs/图3.3.2 分阶段状态区县数分布.png")

